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This paper presents a static and dynamic voltage stability analysis of a power 
network with thermal and wind generation considering static and dynamic 
load models. The thermal plant was modeled as a synchronous machine and 
the wind farm as a variable speed induction generator based on a doubly-fed 
induction generator. The load considered the ZIP, exponential recovery, 
induction motor, and frequency-dependent load models. The bifurcation 
points were found by continuation power flow and sensitivity analyses. In 
addition, dynamic voltage stability assessments were performed considering 
changes in the moment of inertia and the frequency parameters. All 
simulations were carried out in a 4-bus power system and using the power 
system analysis toolbox (PSAT) and MATLAB script code. The results show 
that the thermal generator had difficulties to maintain stability under dynamic 
load variations and frequency changes, the wind generator had difficulties to 


ZIP load model ee ae: i 
maintain voltage for the load with induction motors, and both generators had 
difficulties when the moment of inertia is increased. 
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1. INTRODUCTION 

Power systems are becoming more complex due to the continuous increase in power demand and 
expansion of electrical networks, the inclusion of various types of loads, and the integration of large 
renewable energies and power electronics. Thus, power system operation is subjected to continuous 
disturbance in the generation, transmission, and the load, which causes the reduction of operation margins 
and leads the system close to its stability limits and collapse [1], [2]. Therefore, studying the voltage stability 
(VS) of the power system with different types of generators and loads is required to evaluate the security of 
the network. 

Some researchers have previously studied this problem. In [3], PV curves are used to study the 
integration of wind generation with compensation devices such as SVC in order to determine the maximum 
loadability of the system. In [4], the authors performed a probabilistic voltage stability-constrained optimal 
power flow to improve the voltage stability margin and social welfare. In [5], the authors studied the impact 
of variable-speed wind power generation in long-term VS while operating synchronous generators with 
dynamic reactive power capabilities. In [6], a VS-constrained optimal power flow was presented to determine 
the required loading margin. In [7], static and dynamic reactive power compensation are used to study the 
integration with generators, considering a 10% step disturbance in the dynamic load model. In [8], the 
integration of variable-speed wind turbines on long-term VS is studied; in this case, the capability curves, 
dynamic model of excitation limiters, on load tap changers, and static and dynamic loads were considered. In 
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[9], a doubly-fed induction generator (DFIG) under maximum power-point tracking mode is used to find the 
control effect of load shedding, by sensitivities of L indices in load buses. In [10], a Monte Carlo 
probabilistic methodology was used to evaluate the oscillatory stability margin in a power system with large- 
scale wind generation. In [11], a voltage stability analysis is used to evaluate the effects produced in a power 
network when reactive power is injected in some selected busses. In [12], the performance of a DFIG driven 
by a wind turbine is evaluated using a control vector technique, obtaining that the use of a AC/DC/AC 
converter helped improving stability and reliability of the network. In [13], the optimal location of wind 
turbines in a power grid was considered to evaluate the performance of the network, where the mathematical 
formulation considered objective functions as the maximization of system loadability and the minimization 
of real power loss of transmission lines; besides, some constraints such as security and small signal stability 
were included in the formulation. 

Most studies in the literature show that voltage stability is analyzed by using static load models. 
However, power systems are nonlinear and complex because demand and generation are composed of 
distinct types of elements that generate dynamic behavior, and their combination causes many uncertainties 
in the variation to be studied. There are still no exact models that demonstrate the actual behavior of the loads 
and with constant parameters it is difficult to predict the proximity to the voltage collapse of a power system 
that is subjected to several changes. 

The main purpose of this paper is to study the voltage stability of a power system that considers the 
inclusion of thermal and wind power generation in a network with different static and dynamic load models. 
The aim is to compare the voltage behavior of the system under the maximum load and close to the saddle— 
node bifurcation for the different load models. A sensitivity analysis is performed by increasing the real and 
reactive power of the load and, after determining critical operation close to the voltage collapse, time-domain 
simulations are used to compare the interaction of different types of load with the thermal and wind 
generation. For this purpose, the rest of the paper has been organized as follows: Section 2 describes the 
generator models, load models, and the stability analysis techniques used in the research; Section 3 presents 
the power system used for the test, the results, and the analysis; and Section 4 summarizes the conclusions 
drawn from the study. 


2. METHODOLOGY 

Figure | presents the flowchart of the methodology used for the voltage stability analysis with the 
two types of generators and four types of load models. The procedure consists of integrating to a power grid, 
two generators based on thermal and wind energies, and then the load models, ZIP, ER, IM, and FD, are 
considered. Next, PV curves are calculated by using the continuation power flow for each type of generation 
and each type of load. Finally, time-domain simulations are run to compare the behavior of the power grid 
with the different operation scenarios. 
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Figure 1. Flowchart of the methodology used for the research 


2.1. Generation model 
Most of the electrical generators in power systems are synchronous machines such as hydraulic and 
thermal generators. This is because the rotor uses DC excitation and imposes a synchronous characteristic 
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that allows greater control of the voltage and frequency in the stator. In addition, these types of machines are 
preferred because they use the rotor as field windings to reduce the stress on the brushes as they operate with 
smaller current magnitudes than in generators where the armature is mobile. At the operational level, some of 
the most important variables in this machine are the active and reactive power, moment of inertia, and 
dynamic response parameters such as sub-transient, transient, and steady-state reactances. To evaluate the 
impact of load variations on the generator, these variables limit the voltage and current response at the system 
buses. 


2.1.1. Thermal generator 

Thermal generators are synchronous machines characterized by having a round rotor. The speed is 
greater than that of the salient-pole machine. Therefore, the mass and inertia of these machines are commonly 
greater than other machines, resulting in a longer mechanical start time and being more stable or slow 
forward response against changes produced in the network. Generators are commonly modeled with the d-q 
representation to facilitate the decoupling of the effects between stator and rotor. This allows to model the 
machine related to its windings with an open circuit test and determine the synchronous impedance. Other 
complementary tests such as stand still time response (SSTR) and slip test allow the characterization of the 
transient and mechanical characteristics of the d-q model [14]. 

The thermal generator was modeled as a four order synchronous machine. This type of machine can 
be formulated using the subtransient d-axis voltage ej instead of ej. The differential equation is 
shown in (1): 


Ww Wy; 
n _ (reg +(xq-xq Jig) 
eg = Tl 2 (1) 
qo 
where To is the subtransient time constant of the q-axis in open circuit, xq is the synchronous reactance of 
the q-axis, xq is the subtransient reactance of the q-axis, i, represents the quadrature current, and eg is the 


subtransient voltage of the d-axis. 


2.1.2. Wind generator 

Wind generators use different types of turbines, such as constant speed turbines, variable speed 
turbines (DFIG), and complements between the previous ones [15]. The DFIG facilitates the speed control, it 
is used to damp oscillations of the system caused by the load dynamics in the network or operational 
conditions whose bifurcations generate voltage collapses. It is necessary to verify the wind penetration 
through the dynamic analysis while considering different loading conditions and variations in its parameters 
[16]. The pitch angle control of DFIG is shown in (2): 


a (Kp0(Wm—Wref )—9p) 
6, = ee (2) 
where T, is the pitch control time constant, Ø is the parameter that changed the pitch angle set when W, — 
Wrer exceeds the value +Aw and K, is the pitch control gain. 


2.2. Load model 

There are two major groups for classifying loads: static and dynamic models. The static models 
focus on calculating the updated values of real and reactive powers, which have dependence on voltage and 
frequency values. This category includes the exponential, polynomial, and frequency-dependent models. 
These models must take into account the voltage levels and the changes in the topology of the network that 
affect the characteristics of load behavior. On the other hand, dynamic models are nonlinear such as 
exponential recovery and induction motors models and they can create many variations in the power system 
operation. 


2.2.1. ZIP load model 

The ZIP model can represent real and reactive powers of load as shown in (3) and (4). These 
equations are divided into the represented constant impedance (Z), constant current (I), and constant power 
load (P) [17]. The terms Vo, Py, and Qo are the nominal values of voltage, real power, and reactive power, 
respectively, as determined under normal operating conditions of the power system. Moreover, coefficients a 
and b are used to represent each of type load, where a; and b; are constants used to represent the percentage 
of impedance in the model, a and bz are constants used to represent the percentage of current in the model, 
and a; and b; are constants used to represent the percentage of real and reactive powers, respectively: 
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Qzip = Qo |o, Aj +b, 5 + bs}. (4) 


2.2.2. Exponential recovery (ER) model 

The ER model and its critical effect on voltage stability were studied in [18]. The model represents a 
nonlinear equation that depends on voltage as shown in Equations (5) to (8). Herein, the terms Vj and Po 
represent, respectively, the initial voltage and power before a voltage change, x, is the active power 
recovery, P4 is the power response, T, is the load recovery time constant, at is the transient active load- 
voltage dependence, and as is the steady state active load—voltage dependence: 


TX) = P;(V) — Pa, (5) 

xp = Pa — P,(V), (6) 

P,(V) = P, H. (7) 
at 

PAV) = Po(=) (8) 


2.2.3. Induction motor (IE) 

The induction motor model refers not only to loads where motors are predominating, but represents 
loads whose dynamics during starts or stops require reactive power to sustain its operation. To represent the 
dynamic behavior of the load, the expression shown in Equation (9) is used. The term T; is the mechanical 
torque on shaft [pu], w represents the angular velocity on shaft [pu], and the coefficients of torque are a, b, 
and c: 


T, =at+bw+cw?. (9) 


2.2.4. Frequency-dependent load model (FD) 

The FD model is the general nonlinear exponential voltage-dependent model used in the ZIP load 
[19] as shown in Equations (10) and (11). Herein, the term Aw is the frequency deviation, kp_g is the active 
or reactive power percentage, and « and f are the exponential coefficients of the active and reactive power, 
respectively: 


xP 

=i5G;) Cae (10) 
xQ 

Q = (7) (1 + Aw)2, (11) 


2.3. Continuation power flow (CPF) model 

The CPF model is used in this research to find the PV curve because it is more precise with the 
values obtained in the maximum load points of the system. The methodology uses a predictor—corrector 
vector to find the solution to the modified power flow equations where the load parameter is included [20] 
and, thus, the maximum load point is calculated considering (12) and (13). Herein, the term @is the vector of 
voltage angles, V is a vector of voltage magnitudes, / is a vector that represents the parameters that multiply 
the power of loads, and Ay,, is the maximum parameter that multiplies the power of load in the power 
system. 


F(6,V,A)=0, (12) 


OSASA 


Max" 


(13) 


Thus, real and reactive powers of loads are modified using the parameter 2 as shown in (14), where 
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k is the constant associated with the rate of power change [20]. The term P; is the new real power of load in 
bus i and P? is the initial real power of load connected to bus i. These changes are used to find the maximum 
power that the system can supply to loads and the results can be represented by using PV or PQ curves: 


P =P%(1+Ak). (14) 


In this research, voltage stability analysis is performed by evaluating the behavior of the power 
system after changing some parameters of elements, frequency, or torque. This theory assumes that the 
parameters vary slowly to identify how and when a system becomes unstable [21]. 


3. RESULTS AND ANALYSIS 
3.1. Study case 

Figure 2 shows the 4-bus test feeder used for the test. This network has a main feeder with a voltage 
of 18kV and 50Hz. This power network considers one load connected to bus 3, which is modeled as 
described in the methodology (ZIP, ER, IM, or FD). Two lines are considered in this diagram to separate 
generation and load and they are modeled as a Pl-equivalent circuit. Only one generator will be connected to 
bus 1 at each simulation and the model depends on the type of machine as previously described in 
methodology (wind and thermal generators). A power transformer connects the generator in bus | and it is 
modeled as equivalent impedance. 

With respect to the wind farm, it consists of 25 turbines of 2 MW. The thermal generator has a 
power of 50 MVA, a rating voltage of 18kV, and an armature resistance of 0.01 p.u. as this large power 
machine involves windings of 100mQ. The leakage reactance was neglected, thus the power losses in the 
core were not considered as characteristic of thermal generators associated with its high speed. A higher 
moment of inertia was configured, where the mechanical start time was 12.8s. 

The static ZIP load model was tested with 50% of the impedance factor, 25% of the current factor, 
and 25% of the power factor. The rated active and reactive powers of the load were considered as 
Po = 30 MW and Qo = 15 MVar, respectively. The induction motor in this document has a nominal power of 
20 MVA at 30kV, simulated under a third-order dynamic model. The resistance and reactance of the stator 
are considered as 0.01 p.u. and 0.06 p.u., respectively. The resistance and reactance in the rotor are 
considered as 0.02 p.u. and 0.06 p.u., respectively. The magnetization reactance is considered as 3.5 p.u. and 
the inertia constant is 1.5kWs/kVA. As the IM is modeled as a speed-dependent load that depends on 
mechanical torque, the mechanical torque and the moment of inertia are increased with values between 0.5 
and 5. In the FD, the real or reactive power percentage was modeled by considering « and f with values 
between 0.5 and 10. For all cases, the simulation considers changes in power as 0.5, 0.99, and 1.22 with the 
reference value found in the saddle—node bifurcation. 


Figure 2. System model of the test that considers the generation and loads models previously described [17] 


3.2. Steady state analysis 

The saddle—node bifurcation was identified by using the continuation power flow method to find the 
maximum load margin around the nominal power load, considering the different load types. The results are 
shown in Table 1. 
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Table 1. Load Margin Considering Load Type 


Load type Load margin 
ZIP 2.555 
ER 2.706 
IM 2.706 
FD e 


After identifying the saddle point as the voltage stability limit, some power variations are performed 
at 0.5, 0.99, and 1.2 times the maximum power, with the aim of identifying the dynamic behavior of voltage 
magnitudes in the network. These simulations are carried out by using the CPF to obtain the curves shown in 
Figure 3. This figure considers the load model connected to bus 3, the power network connected to bus 4, and 
PV generation connected to bus 1. 

Figure 3(a) shows the PV curves with a ZIP model. The results show that an increase in the system 
loadability leads to a reduction of the maximum load parameter À calculated with the CPF. Thus, a smaller 
power variation leads the power system to a fast voltage collapse because the system is operating close to the 
voltage limits. 

Figure 3(b) shows the PV curves with the ER model or the IM model. The results show the PV 
curves are similar for both load models, as the continuation power flow does not consider dynamic 
parameters of load and the increase is performed with static steps. The simulations with the FD model were 
not succeeded because this model is a frequency-dependent load and the steps of the algorithm do not 
consider these variations. 
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Figure 3. PV curves for different power loads for the (a) ZIP model and (b) ER and IM models 


3.3. Time-domain simulation 

In [17], some simulations were performed to identify several Hopf bifurcations by using the same 
test case. In this paper, we focus on the interpretation of the different behaviors of the voltage magnitude with 
the time-domain simulation, starting from the critical values of the load in the power system. Thus, Figure 4 
shows the behavior of the voltage magnitude in buses 1 and 3 when the nominal load is multiplied by 
1.8 times and two types of load are connected to bus 3. 

Figure 4(a) shows the response of the system under this operation condition when the load is a ZIP 
model. The results show that both buses begin with a high voltage magnitude but when the load increases, the 
DFIG overspeeds and the voltage drops in about 2 seconds. Although the voltage dropped, it did not collapse 
and the system is capable to stabilize. 

Figure 4(b) shows the voltage collapse with the IM load model when the real power is multiplied by 
the same load factor, as previously defined. This figure shows the bifurcation points when the power changes 
with a lower value than the maximum power determined previously with the CPF. This point can be 
associated with the Hopf bifurcation point, because the voltage magnitudes of the buses change abruptly and 
then reach zero in a few milliseconds. 
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Figure 4. Voltage magnitudes calculated with time-domain simulations considering (a) ZIP and (b) IM 


Figure 5. Time-domain simulation considering different power values with (a) ZIP model and wind 
generator, (b) ZIP model and thermal generator, (c) ER model and wind generator, and 


Figure 5 shows the voltage magnitudes in bus 3 calculated with time-domain simulations. In these 
cases, voltage magnitudes are compared when the ZIP or ER load models are connected to the bus 3 and the 
thermal generator or the wind generators are connected to bus 1. Figure 5(a) shows a voltage collapse when 
the real power is multiplied by a factor of 1.2, whereas Figure 5(b) shows that the system does not reach 
voltage collapse. For the thermal generation case, the voltage magnitudes are higher than those obtained with 
the wind generation case. Figure 5(c) and Figure 5(d) show that the ER load does not present as inconvenient 
during the operation for the considered thermal and wind generator, i.e., both are capable of maintaining 
voltage stability despite of the factor multiplied to reach the power limit. 
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(d) ER model and thermal generator 
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Figure 6 shows the voltage magnitudes in bus 3 calculated with time-domain simulations. In these 
cases, the IM load is connected to bus 3 and the inertia moment and the torque are changed. The results show 
that voltage collapses are presented in the system when wind generators are considered in the study. 
Figure 6(a) and Figure 6(b) show that the voltage magnitudes have different behaviors considering the torque 
variation. Figure 6(a) shows that the system collapses when the torque increases by a factor of 0.99 and 
Figure 6(b) shows that the voltages in the three variations have the same behavior and the system stabilizes. 
Figure 6(c) and Figure 6(d) show large voltage magnitude variations in bus 3, meaning that thermal and wind 
generators are more sensitive to the increase of inertia. It is noted that the wind generator is unable to respond 
to a factor that increases seven times, resulting in a voltage collapse; however, the thermal generator 
oscillates around 10 seconds to the same event, but stabilizes after few seconds. 
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Figure 6. Time-domain simulation considering IM load model with (a) wind generator and torque variation, 
(b) thermal generator and torque variation, (c) wind generator and inertia moment variation, and 
(d) ER model and inertia moment variation 


Figure 7 shows the voltage magnitudes in bus 3 calculated with time-domain simulations. In these 
cases, the FD load is connected to bus 3 and the real or reactive power percentages are changed. The results 
show that the system is very sensitive to parameter variations and both generators are not able to respond to 
the load percentage with the FD model and system collapses. 

The results show that there is no difference in the PV curves obtained by the CPF for wind and 
thermal generators when dynamic models are considered, because the CPF does not reflect the dynamic 
behavior for the static characterization of the PV curve. When nominal power increases closer to the saddle— 
node bifurcation, the changes around the nominal power are smaller and a reduction of the load factor is 
presented. Thus, Table 2 summarizes all the results obtained with the time-domain simulations with different 
parameter variations. 
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Figure 7. Time-domain simulation considering an FD model for (a) wind generator and (b) thermal generator 


Table 2. Summary of the Results Obtained from the Time-Domain Simulations 


Load type Parameter Wind generator Thermal generator 
Nominal Voltage reduction as it approaches the power limit Although a voltage reduction is presented, the 
ZIP ous and voltage collapse with a greater power increase thermal generator is more capable to maintain 
P after the maximum loadability. stability. 
ER Nominal Voltage magnitudes increase on the bus, but the No instabilities are presented for the considered 
power system remains stable. cases. 
ree ie : Volt scillations stabilized ickl 
; Oscillations are stabilized without overvoltage but E S ara Hee Neate 
Mechanical À ; : but overvoltages and reductions of the final voltage 
IM may result in collapses or reductions in the voltage A 
Torque ; level are generated. However, under simulated 
level in the bus. 
cases no voltage collapse occurs. 
Overvoltages are generated during the voltage 
M Moment of A voltage collapse occurs when the initial moment oscillation, being stable after 15 or 20 seconds for 
inertia of inertia is increased. the simulated conditions. There is no voltage 
collapse. 
Thi iations i d fr t F : 
OO rans erat nen ee The voltage dynamics are not sustained by the 
generate oscillations and voltage collapses themial ocneration. when the load depends ‘on 
F Kp quickly. Therefore, it is insufficient to handle the g pone 


inertia of the generator to compensate for the 


frequency changes, even being a generator with a 
greater moment of inertia. 


variances of frequency demanded by this load. 


4. CONCLUSION 

A static and dynamic stability analysis was performed in a 4-bus system with different generation 
and load types to identify the voltage behavior with parameter variations. PV curves were calculated to find 
the maximum loadability of the power system (saddle point). Additionally, we performed time-domain 
simulations considering wind and thermal generators and different load models because PV curves do not 
consider the type of generation as it is treated as a PV bus. From the results, we found that the exponential 
recovery load has the widest range of variation in the different cases and improves the voltage magnitudes. 
With the frequency-dependent load, the margins were unknown because of dynamic load models are not 
considered by the CPF model. A voltage collapse occurred with the ZIP load model when the power was 
increased beyond the limit because this static load model allows creating the PV curves. The induction motor 
generates collapses when increasing the moment of inertia. The ZIP load model with thermal generation does 
conduct the system to a voltage collapse with the simulated cases and the ER model does not present 
instabilities with both generators. The motor has a better response to voltage fluctuations and the system is 
capable to maintain stability with a thermal generator. Finally, the FD load model generates many 
singularities in the solution at the moment of calculating the response with the time-domain simulation. 
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